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The  purpose  of  this  Investigation  is  to  find  the  theoretical 
solutions  to  scattering  of  electromagnetic  fields  by  dielectric  bodies 
of  revolution.  This  is  the  first  part  of  the  final  goal  of  solving 
scattering  by  burled  dielectric  objects.  The  principal  technique 
used  in  this  investigation  is  the  unimoment  method  which  was  developed 
for  the  purpose  of  solving  scattering  by  dielectric  bodies  under 
contract  DAAK02-71-C-0206.  The  application  of  this  method  to  scattering 
by  dielectric  cylinders  in  a 2-dimenslonal  space  has  been  reported  in 
the  final  report  of  the  above  contract.  In  this  report,  we  shall  extend 
the  method  to  scattering  in  3-dimensions.  In  the  application  of  the 
unimoment  method,  one  frequently  has  to  deal  with  the  finite  difference 
or  finite  element  method.  In  this  report,  we  start  with  a discussion 
on  these  two  methods  regarding  their  differences  and  similarities.  A 
method  of  Inverting  banded  sparse  matrix  Is  also  developed,  and  finally, 
the  development  of  a special  potential  formulation  of  the  Hixwell's 
equations  for  the  axial  synmetrlc  problems,  and  the  solution  of  the 
scattering  problems. 
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ABSTRACT 

Recent  development  in  the  unimoment  method^  has  brought  finite  difference 
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and  finite  element  methods  into  the  computational  techniques  of  electro- 
magnetic scattering.  In  this  paper  the  various  finite  methods  and  their 
potentials  in  the  scattering  computations  are  examined.  A section  on 
programming  technique  is  Included  for  those  uninitiated  and  the  applications 
of  the  finite  methods  in  2-dimensional  and  3-dimensional  scattering  problems, 
together  with  some  of  tne  associated  computational  subject  matters  are 
presented . 
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Introduction 


Recent  development  In  the  unimoment  method  has  brought  finite 
2 3 

difference  and  finite  element  methods  into  the  computational  techniques 
of  electromagnetic  scattering.  Finite  difference  had  been  a dominant 
computational  method  of  solving  partial  different!  1 equations,  until  it 
was  overshadowed  by  the  finite  element  methods  in  the  mid-sixties.  Recent 
development  in  sparse  matrix  algorithms  has  also  changed  the  manner  by  which 
the  finite  methods  are  Implemented.  Originally,  the  relaxation  methods  were 
the  main  tools  of  solving  finite  difference  equations,  but  now  almost  all 
finite  equations  are  Inverted  Instead  of  iteratively  solved.  The  finite 
element  method  which  was  first  developed  by  the  struc tura 1 engineers  is  now 

applied  in  almost  all  engineering  disciplines.  It  is  now  well  understood. 
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Mathematical  analysis  and  sophisticated  techniques  ’ have  been  developed 
for  it.  In  antenna  and  scattering  computations  the  finite  techniques  are 
still  like  new  found  friends,  and  the  acquaintance  with  them  is  yet  to 
be  developed.  The  purpose  of  this  paper  is  to  examine  those  methods  from 
the  stand  point  of  one  who  is  interested  in  electromagnetic  scattering 
computations.  A section  on  programming  technique  is  included  for  those 
uninitiated  and  the  applications  of  the  finite  methods  in  2-dlmensional 
and  3-dlmenslonal  scattering  problems,  together  with  some  of  the  associated 
subject  matters  are  presented. 

Finite  Methods 

In  this  paper  the  term  finite  methods  is  meant  to  include  finite 
difference,  weighting  function,  Galerkin's  and  finite  element  (or  Ritz 
finite  element)  methods.  Literally  any  numerical  method  that  approximates 
continuffl  mathematics  by  discrete  mathematics  may  be  termed  a finite  method, 
however,  that  is  not  usually  what  we  mean.  For  Instance  the  moment  method 
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of  solving  integral  equations  is  not  classified  as  a finite  method.  A 
finite  method  is  traditionally  used  to  solve  differential  or  partial 
differential  equations  directly. 

During  the  last  fifteen  years  computation  in  electromagnetic  scattering 
problems  have  been  actively  pursued  almost  entirely  in  terms  of  integral 
equations.  Needless  to  say  the  drift  to  integral  equation  is  very  natural 
and  for  good  reasons.  In  integral  equations  the  computations  are  limited 
to  the  scatterer  or  antenna  while  in  tlie  finite  methods  they  are  generally 
spread  over  the  entire  space.  In  integral  equations  the  radiation 
conditions  are  automatically  satisfied  while  in  the  finite  methods  they 
require  special  numerical  treatments  which  are  often  unsatisfactory.  Recent 
advances  in  communications  and  remote  sensing  have  urgently  demanded  the 
results  of  scattering  by  dielectric  and  lossy  inhomogeneous  bodies.  The  only 
practical  approach  to  such  problems  appears  to  be  direct  solution  of 
partial  differential  equations  rather  than  solution  by  integral  equation, 
the  formulation  of  which  in  an  inhomogeneous  medium  is  already  i difficult 
task.  Recent  development  in  the  applications  of  finite  methods^  has 
essentially  cleared  the  way  for  general  applications  of  the  finite  methods 
in  electromagnetic  scattering. 

For  the  sake  of  brevity  we  shall  discuss  the  methods  first  in  one 

dimension.  Let^be  a differential  operator  or  more  specifically  for 

2 

Illustration  we  set  ^ . The  objective  of  the  finite  methods 

dx 

is  to  solve  the  equation 
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In  most  numerical  methods,  the  independent  variable  is  discretized 
as  shown  In  Fig.  1.  The  formulations  of  the  finite  methods  are  illustrated 
in  the  following: 


Finite  difference  from  fundamental  meaning  of  derivative: 


^+1  ■ ^-1 


= f . i = 1 , . . . ,N-2 
1 


Weighting  function  method  I (Wra  I) 


f W.(x)dx  = f f(x)W  ( 

Jo  dx^  ^ Jo  ^ 


x)dx  i=l,2...M 


Weighting  function  method  II  (WFM  II)  (Integration  by  parts  of  WFM  I) 


A I dW,  1 

<Ji'(x)W,(x)  - f dx  = f f(x)W.(x)dx 

^ 0 Jo  Jo 


Galerkln's  method 


4'(x)  = ^ a.t.(x)  and  set  W,(x)  in  (4)  to  ;,(x) 

i = 1 1 ^ 1 

(Ritz)  Finite  Element 

Minimize  I = f [ (t^)  ^ + 2fiJ)]dx 

Jo 


In  most  text  books  the  weighting  function  methods  are  not  distinguished  by 
whether  integration  by  parts  is  used  or  not.  For  convenience  of  the 
following  discussion  we  shall  c.all  it  weighting  function  method  II  when 
integration  by  parts  is  used,  otherwise  weighting  function  method  I. 

We  shall  name  the  set  of  functions  which  are  used  to  approximate 


the  solution  li(x)  the  trial  functions,  thus 


(6) 


N 

4i(x)  = V a r . (x) 
i = 1 

where  are  trial  functions.  W^(x)  are  a set  of  weighting  functions. 

The  following  are  a few  definitions  of  commonly  used  terms: 

patch  function  - functions  which  are  zero  outside  an  element  and  its 
close  neighbors,  such  as  those  shown  in  Fig.  2. 

element  - a line  segment  (one  dimension),  area  (two  dimension)  or 

volume  (three  dimension)  as  a result  of  discretization. 

interpolate  - a set  of  trial  functions  is  termed  interpolat i ve 

if  it  is  unity  at  one  node  and  zero  at  all  other  nodes.  An  interpolat ive 

trial  function  has  the  property  Chat 
N 

♦(x)  ' V a ; (x) 
i = 1 

where  represents  fi(Xj). 

conformity  - a set  of  trial  functions  is  termed  conforming  if  each 
of  its  members  Is  continuous. 


Only  those  methods  that  use  patch  functions  as  trial  functions  and 
weighting  functions  will  be  considered  as  finite  methods. 

It  is  well-known  that  the  finite  difference  equations  are  derived  from 
Taylor's  series.  It  is  not  difficult  to  relate  the  finite  difference 
equation  to  the  W.F.  method.  Consider  a quadratic  patch  function 


Q^(x)  = + B^(x  - x^)  + C^(x  - Xj^)‘ 


passing  through  <t'(x^_^),  il>(x^),  ij>(x^^^).  It  is 
\ :C.= 


2Ax 


i 


(7) 

immediately  evident 

- 2}!(x^)  + 't'(Xj^_j) 

2(Ax)^ 
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In  (3)  let  4>(x)  = ^ Q.  (x)  , tlie  weighting  tiinctions  W.(x)  be 

i = L 

6{x  - X.),  and  we  tliid  that  (3)  lediu  es  to  the  fltilte  difference  equation  (2). 

^ M 

On  the  other  hand,  in  (5)  let  i>(x)  = y'  a.f,  .(x) 

1 1 

1 = 1 


-i  M '7 

1 = 1 U (x)  ] + 2f  53  a ' (x;}dx 

Jo  i = i i = 1 

1 M 


M. 

3a 


X.  ' (x)  If,  ' (x)  + 2ff,.  (x)  )dx  = 0 
■ ' J 1 


j = 1,2, ...M 


(H) 


which  is  equivalent  to  (4)  if  the  trial  functions  f^(x)  and  the  weighting 
funcf'  s W^(x)  are  the  same  set  of  functions,  and  W^fO)  = “ O- 

r ^_t  r i cjii  o n^ 

While  W.F.  method  I and  II  seem  to  be  univ  a simple  mathematical 
step  away  from  each  other,  their  numerical  constraints  are  quite  different. 
In  the  simple  example  of  Eqs.  (3)  and  (4),  It  is  evident  that  in  W.F. 
method  I the  trial  functions  must  be  at  least  Quadratic  and  the  weighting 
functions  be  integrable;  while  in  W.F.  method  II,  the  trial  functions  can 
be  linear  and  the  weighting  function  must  belong  to  C',  (W (x)  integrable). 
The  consequence  is:  one  can  use  linear  functions  as  trial  functions  in 

W.F.  method  II,  but  not  in  method  I. 

Another  important  observation  is  that  the  union  of  the  set  of  trial 
functions  must  cover  the  entire  region,  while  the  weighting  functions  are 
not  so  required.  Thus,  in  the  W.F.  methods,  the  weighting  functions  may 
be  sampling  functions  e.g.  pulses  or  patch  functions  which  do  not  overlap, 
such  as  shovm  in  Fig,  3 as  oppose*!  ti'  those  patch  functions  of  Fig.  2 


which  are  used  as  trial  functions. 


FT 


1 

One  important  property  of  sampling  is  that  discrepancies  in  the  trial 
functions  may  be  suppressed  by  the  weighting  functions.  An  obvious 
discrepancy  in  the  W.F.  method  1 using  piecewise  quadratic  trial  functions 

is  that  the  trial  functions  are  not  conforming,  as  illustrated  in  Fig.  4. 

Ax  Ax 

If  we  use  the  left  of  the  quadratic  for  x.  <x<x.+—  and 

1 2 — — 1 2 

the  right  quadratic  for  - ~ .1  x ^ ^i+1  ^ discontinuity 

is  obvious.  Yet  It  causes  no  ill  effect  because  the  weighting  functions 
are  zero  at  the  discon' inuities. 

The  case  of  discontinuity  becomes  very  serious  for  the  finite  element 
method,  since  the  trial  functions  are  now  chosen  to  be  the  weighting 
functions  and  the  union  i>i  trial  i uni  Lions  must  cov.  r the  entire  region.  The 
derivatives  at  the  discontinuities  result  in  delta  functions,  and  the 
square  of  delta  functions  are  not  integrable.  That  is  why  elements  must 
be  conforming  in  the  finite  element  method.  However,  being  non-conforming 
may  still  not  be  fatal  in  the  finite  element  method.  Frequently  the 
irregularities  are  ignored  and  the  computations  still  may  yieli  good 
results,  although  just  as  frequently  they  may  fail.  Nonconforming  elements 
are  often  referred  to  as  "variational  crime"  in  finite  element  literatures, 
and  there  exists  a test,  known  as  the  "patch  test,"  which  may  be  used 
to  predetermine  whether  a specific  sot  of  non-conf orming  elements  would 
converge  to  the  correct  result.  In  one  dimension,  conforming  elements 
are  easily  obta ined  even  for  high  order  polynomial  'approximations.  Any 
piecewise  polynomial  function,  which  is  non-zero  at  only  one  node  is 
a conforming  element.  In  two  or  more  dimensions  conforming  requires 
extra  effort,  which  may  result  in  increasing  the  number  of  nodes  or  the 
hand  of  the  matrix  or  both. 
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One  Dimension  Fin i M^t tic^ s 

Consider  the  W.F.  method  II  using  linear  functions  for  both 
C^(x)  and  W^(x)  such  as  those  shown  in  Fig.  2a.  The  integral  in  (4) 
from  becomes 


Ax 


- 2fi.  + 


i-1 


":X 


Ax 


Ax 


dx 


(9) 


which  is  almost  identical  to  the  finite  dilference  equation  (2)  except 
that  the  right  hand  side  of  (2)  is  collocation  and  that  of  (4)  is 
weighted  by  a linear  patch  function.  Indeed  the  base  of  the  weighting 
function  may  be  reduced  without  any  effect  on  the  equation. 

Now,  consider  a quadratic  trial  tunction  used  in  W.F.  method  II,  i.e.. 


q^(x)  = 


--  (x  - X.)  + ^ (x  - X.) 

^ 2 Ax"  ^ 


(10) 


then. 


d?.(x) 

dx 


^i+1  ■ ^i  , ^+1  ‘ 2^i  ^ n-1  . 
"1^'  ^ — :t - >^i) 

Ax 


(11) 


dW. 

If  the  weighting  function  is  even  with  respect  to  x^,  will  be  odd 

function  w.r.t.  x.,  thus  the  first  term  in  (11)  will  have  no  contribution 
^ dC^(x)  dW^(x) 

in  the  Integra)  of  the  product  — ^ , and  we  shall  get  an 

equation  again  identical  to  the  finite  difterence  equation. 

It  seems  that  all  finite  methods  can  be  reduced  to  the  same  equations. 
The  fact  l.s  that.  In  this  particular  example,  the  second  order  polynomial 
operated  on  by  tlu'  differential  I'pcrator  results  in  a constant.  As  a result. 


i 
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for  this  example  a sampling  weighting  function  set  Is  as  good  as  any 
distributed  weighting  function  set.  Indeed,  the  quadratic  functions  are 
the  best  approx lm.at  Ion  functions  one  can  get  to  solve  this  partliular 
problem,  because  the  differential  equation  Is  satisfied  not  only  at  the 
nodes  but  also  between  nodes,  and  any  attempt  to  use  higher  ordered 
approximations  usually  makes  the  approxinvit ion  worse  rather  than  better. 

Two  Dimensional  Problems 

The  above  is  the  story  for  one-dimensional  equations.  In  two 
dimensions  many  of  the  above  discussions  are  still  valid,  but  there  are 
other  complications. 

(A)  Conformity 

Conformity  of  elements  In  one  dimension  is  never  a problem,  *'jt 
in  two-dimensions  it  is  a major  problem  when  liigh  order  approximations  are 
to  be  considered.  Consider  a typical  finite  element  discretization  in 
Fig.  5.  The  elements  are  triangular  and  the  trial  functions  are  pyramidal 
like  linear  patch  functions.  Obviously  a visual  representation  of  the 
relation  between  our  trial  functions  and  the  problem  geometry  requires  the 
use  of  one  addltlon.al  dimension.  A Lvplcal  trial  I unction  consists  of  6 elements 

and  7 nodes.  Referring  to  Fig.  5,  a linear  trial  function  centering  at 
the  point  i,  consists  of  six  planes  which  have  values  one  at  the  point  1 
and  0 at  the  other  nodes,  such  as  shovm  in  Fig.  b. 

It  is  evident  that  the  linear  elements  are  conforming,  i.e.,  any 
linear  combination  of  such  trial  functions  shall  have  no  discontinuity. 

The  second  order  approximation  in  the  finite  element  methods  uses 
a piecewise  quadratic  function  for  t(x,y)  within  each  element, 


2 2 

)(x,v)  = a„  -t-  a,x  + a_y  + xv  + a.x  + a,v 
0 121  a 5' 

-q- 


a2) 


The  quadratic  function  of  (x,y)  has  six  dcwrees  of  freedom.  The 
three  nodes  of  a triangular  element  are  not  sufficient  to  determine  the 
six  coefficients  of  the  quadratic.  Extra  nodes  need  to  be  added,  such 
as  shovm  in  Fig.  7,  in  which  a node  is  designated  on  each  side  of  an 
element.  These  trial  functions  are  now  conforming.  (Quadratics  which  are 
unity  at  one  node  and  ^ero  at  the  rest  of  the  nodes  of  a triangular 
element  can  be  found  either  by  inverting  a 6 x 6 matrix  or  by  using  an 
area  coordinate  system  (for  details  see  refs.  (5),  (7j.)  The  second  order 
approximation  in  the  finite  element  method  gives  better  accuracy  but  it 
also  increases  the  number  of  nodes  and  increases  the  bandwidth  of  the 
matrix. 

In  finite  difference  assuming  square  mesh,  the  usual  five  point 
equation  for  the  Laplacian  operator  gives,  (referring  to  Fig.  5) 


It  is  a result  of  quadratic  approximation  of  (f(x,y)  near  but  it  is 

not  a complete  quadratic.  Indeed,  a five  point  equation  cannot  represent 
a complete  quadratic  which  has  six  coefficients.  The  missing  term  in 
(13)  is  a^  X y,  which  can  be  omitted  if  the  mesh  is  rectangular. 

The  trial  functions  result  in  (13)  are  not  conforming  but  the 
discrepancy  is  again  suppressed  by  the  weighting  functions.  Again,  for 
solving  Laplace's  equation, the  quadratic  equation  of  (13)  satisfies 
the  differential  equation  between  nodes  along  the  mesh,  and  it  should  be 
a better  approximation  than  the  nine  point  finite  difference  equations 
(fourth  order  polynomial)  which  satisfies  the  differential  equation  only 
at  the  nodes.  Many  engineers  have  been  surprised  to  find  that  their 
nine  point  finite  difference  equations  give  worse  results  than  the 
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five  point  finite  difference  equations  when  dealing  with  Laplace’s  equation. 


(B)  Consistency 

Because  of  the  sampling  property  of  the  finite  difference  equations, 
there  are  frequent  inconsistancies,  particularly  when  dealing  witii 
inhomogeneous  media.  In  the  case  of  a Helmholtz  equation  for  a 
continuously  inhomogeneous  medium,  the  finite  difference  equation  is, 

^+N  ^+1  ^ ^-1  ® 

where  the  field  points  ■ . <t>.  0.  ,,  are  assumed  to  be 

i+N  1+1  1 i-1  i-N 

in  a homogeneous  medium  Let  1 be  replaced  by  i + 1 , we  find 

are  now  assumed  to  be  in  a homogeneous  medium  This  type  of 

inconsistency,  while  tolerable,  is  not  very  satisfying  theoretically. 

The  problem  becomes  more  serious  when  a node  is  situated  at  a material 
discontinuity  such  as  shown  in  Fig.  8.  There  are  several  ways  of 
handling  material  discontinuities  in  the  finite  difference  method  but 
in  general  they  are  not  satisfying. 

In  the  finite  element  method,  the  fields  are  clearly  identified  in 
each  patch,  hence  there  is  no  ambiguity  like  that  of  the  finite 
difference.  Furthermore,  the  quadratic  functions  in  the  finite  element  are 
piecewise  quadratic  and  a typical  trial  function  actually  consists  of 
six  piecewise  quadratic  functions,  while  in  the  finite  difference  method 
the  requirement  of  twice  differentiability  on  the  trial  functions 
necessitates  a single  polynomial  to  cover  a group  of  connected  nodes. 

(C)  Flexibility 

Because  the  polynomials  in  the  finite  elements  method  are  complete, 
the  method  is  essentially  indepeiuieiU  of  the  shape  of  the  triangular 
element.  In  ttic  llniio  d i f feriMu  i’  nn-tiiod,  tlic  polvnomlal  is  Incomplete, 


i 


hence  it  is  restricted  to  rectangular  meshes.  It  can  be  made  to  be 
independent  of  mesh  geometry,  however,  bv  using  a six  point  difference 
mesh  such  as  showm  in  Fig.  9.  Indeed,  for  a six  point  finite  difference, 
since  it  is  now  independent  ot  mesh  geometry,  tlie  mesh  is  no  longer 
needed.  The  only  requirement  is  that  the  nodes  are  relatively  uniformly 
distributed.  A survey  of  recent  developments  in  variable  grid  finite 
difference  methods  can  be  found  in  Ref.  (8). 

The  usual  applications  of  ttie  finite  difference  method  use  rectangular 

mesh,  with  the  boundary  nodes  specially  treated  bv  extrapolation  or 
9 

interpolation.  The  merit  of  such  a shceiae  is  the  simplified  programming, 
and  the  disadvantage  is  reduced  accuracy  and  the  difficulties  involved 
in  finding  the  normal  derivatives  from  the  solution.  In  the  finite  element 
method,  using  triangular  elements,  it  is  necessary  to  store  information 
about  the  locations  of  the  nodes,  and  also  the  information  about  the  elements, 
i.e.,  the  nodes  which  form  <t  particular  element.  Finite  difference 
raethodsusing  variable  mesh  should  require  the  stime  information,  i.e., 
the  information  about  the  neighboring  nodes  ot  a control  node  must  be 
available.  Consequently,  finite  difference  method  using  variable 
mesh  is  about  as  flexible  as  the  finite  element  method,  and  the  degree 
of  complexity  in  programming  for  both  methods  are  about  the  same. 

The  conclusion  is  that  the  finite  difference  method  is  about  as  good 
as  linear  approximation  In  the  finite  element  method.  The  higher  order 
approximation  in  the  finite  element  method  should  be  more  accurate  than 
the  corresponding  approximations  in  the  finite  difference  method,  yet  the 
finite  element  method  would  require  more  nodes.  The  finite  difference 
method  using  a rectangular  grid  Is  simple  to  program  but  less  flexible,  and  by 
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using  star  meshes  it  can  be  made  more  flexible  but  its  complexity  is 
raised  to  that  of  the  linear  finite  element  method.  The  method  not  very 
often  used  is  the  W.F.  method  II,  with  which  one  may  use  piecewise 
polynomial  trial  functions  as  in  the  finite  element  method  yet  with  a 
much  relaxed  conformity  requirement.  It  will  be  a worthwhile  project  to 
investigate  whether  the  relaxed  conformity  requirement  on  W.F.  method  II 
can  reap  some  of  the  benefits  of  higher  order  approximations  without 
increasing  the  nodes  or  bandwidth  of  the  resulting  matrix. 

Programming  Technique 

To  implement  the  finite  methods  one  need  first  to  describe  the  geometry 
and  construct  the  mesh,  followed  by  generating  the  matrix  and  solving  the 
matrix.  These  basic  steps  are  described  in  the  following; 

( 1 ) Mesh  Const ruction 

In  most  practical  problems,  the  node  numbers  are  large  and  the 
mesh  topology  is  complicated.  Therefore,  it  is  almost  an  impossible  task 
to  instruct  the  computer  In  an  exact  mesh  geometry.  The  best 
strategy  is  to  construct  a rule  which  the  program  can  follow  to  generate 
the  mesh  au toraa t ical ly . There  are  many  ways  of  automatic  mesh  generation. 

In  the  following  we  shall  describe  one  of  them,  which  we  have  used 
very  effectively: 

( a ) Construct  a regular  mesh  ins ide  t he _j^l r c ular  boundary 

This  can  be  done  in  many  ways  a few  of  whicli  are  shown  in  Fig.  10. 

These  mesh  can  be  generated  automatically.  Fig.  10(a)  is  used  if  the 
center  part  of  the  circle  is  Included.  Fig. 10(c)  is  used  when  the  center 


part  of  the  circle  is  not  Included  in  t)ie  finite  method. 
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(b)  Find  the  iiiters^ctiojis  of  scatterer  contour  with  the  mesh 

line 

If  the  mathematical  formula  of  the  scatterer  contour  is  known,  it 
is  relatively  easy  to  find  the  intersections  of  the  contour  with  tlie 
mesh  lines.  For  the  mesh  geometry  of  Fig.  10(c)  it  is  a particularly 
easy  task,  because  we  only  need  to  find  the  intersections  between  the 
contour  and  the  radial  lines.  One  of  the  two  nearest  nodes  to  an 
intersection  point  along  the  radial  line  is  now  shifted  to  the  intersection 
point.  The  rest  of  the  nodal  points  along  the  particular  radial  lines 
are  now  shifted  so  as  to  make  the  perturbation  equally  distributed.  It 
is  noticed  that  there  are  two  possible  orientations  of  the  diagonal 
hypotenuse  of  the  triangular  elements  between  two  adjacent  radial  lines. 

We  shall  call  those  elements  of  Fig.  11a,  left  handed  (LH)  elements  and 
those  of  Fig.  lib,  right  handed  (RH)  elements.  Ihe  choice  on  the  type 
of  element  may  be  made  based  on  the  rule, 

R, . , > R.  =*  LH  element 
i+1  — 1 

R, ,,  < R.  ^ RH  element 
i+1  1 

Using  the  above  rultis,  the  contour  of  a distorted  raindrop  is  discretized 
in  Fig.  12. 

(2)  Computation  of  Matrix  Elements 

The  key  to  the  programming  of  any  finite  method  with  variable  elements 
is  the  name  list  of  nodes,  which  form  a particular  element.  Referring  to 
Fig.  13,  we  generate  an  array  ELKM(I,J),  where  1,  Is  the  numbering  of  the 
element  and  the  dimension  for  J is  three.  This  array  is  used  to  store 
the  list  of  the  nodes  which  form  the  particular  element  1.  For  example, 
F.LEM(12,1)  = 3;  ELF.M(12,2)  = 6,  ELf><(12,3)  = 7.  Of  course,  the  coordinate 
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locations  of  each  node  should  be  stored  in  another  array. 


With  the  information  of  the  elements  and  nodes  available,  we  now 
consider  the  generation  of  the  nuitrix  elements.  In  the  following,  we 
shall  consider  the  generation  of  the  matrix  element  of  the  two-dimension.al 
Helmholtz  equation  using  the  weighting  function  method  II.  Thus  we 
compute  the  integral. 


(VU  • 7W.  + k UW.)ds 
11 


Both  U and  W will  be  piecewise  linear  functions,  therefore  the  result 
will  be  identical  to  the  finite  element  method. 

A typical  linear  element  h,  with  nodes  (h)j^,  (h)^,  (h)^  should  have 
the  field  within  the  element  described  by  a linear  function. 


^h)  = ^h)i  ^h)2"  ^ ^h)3y 


The  coefficients  a are  found  by  matrix  inversion, 

th)  1 


^(h)  ■ ^(h)  ^(h) 


The  vector  L'  , , consists  of  the  field  at  the  nodal  points  hence  U.,  .... 

(h)  in;i 


(i  = 1,2,3) , and 


^ ’‘(h)l  ^(h)l 

^ ^(h)2  ^(h)2 

^ ’‘(h)3  ^(h)3 


Let  us  consider  W.  to  be  a linear  pyramid  function  which  is  unity  at  node  1 
1 

and  zero  at  tlic  surrounding  six  nodes,  such  as  shown  in  Fig.  14.  Throughout 
the  entire  process  of  weighting  integrals,  each  element  will  be  weighted 
by  three  different  weighting  functions,  each  of  which  may  be  represented 
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by  a linear  function  which  is  unity  at  one  node  of  the  element  and  zero 

at  the  other  two.  The  function  whicli  is  unity  at  node  (h)i  is  actually 

part  of  the  weighting  function  (ii)i,  and  the  coefficients  generated  by  it 

are  associated  with  the  (h)i^^  equation  or  the  (h)  i row  of  ttie  rnatrLx. 

To  each  in  an  element  h,  there  are  three  associated  weighting 

(h) 

functionsW^^^.  (j  = 1,2.3).  Let  be  part  of  which  is  associated 

with  the  element  h and  whicli  has  value  at  the  nod'^  (h)j  and  zero  at  the 
other  two  nodes  of  element  h. 


(h)l  ‘'(h)  2 


Hm  ■ "j 


where 


- N - "1  - [° 

n = PI  , n = 1 , n = 0 

^ loj  lo  llJ 


The  gradient  of  gives; 


0 0 0 

= 010 
0 0 1 


0 0 0 
0 1 0 
0 0 1 


0 I 0 N 1' 

on,  (h)  (h) 


and  the  gradient  of  ”(h) 

1 P 0 _ 

TW^..  = 0 1 0 N n. 

(h)  [o  0 J (h)  J 

The  Integral  | '^^'(h)  * '^(h 


VW-;,  , ds  is. 


0 

0 

0 

-1  T 

0 

0 

o' 

0 

1 

0 

0 

1 

0 

0 

0 

1 

(h) 

0 

0 

1 

where  is  the  area  of  the  element  h,  and  P(h)3^  ' 

A . xP"!  is  recognized  as  the  corn  r ibnt  ion  to  the  matrix  element 
(h)  (h)i 

G[(h)j,  (h)i]  (i  = 1,2,J).  This  value  should  be  accumulated  in  the  array 
G in  the  location  C((h)j,  (h)!). 

Simularly  the  second  integral  of  (15)  can  be  carried  out  by 


I , ^h)  "(h)  ^h)  " 


(1  X v)  N,u\  n.ds 
(h)  ] 


1 X V 

r.'l  T 2 .,-1  - a 

(h)  (h)  2 (h)  ] 

V XV  X . 


2 r r.T 


<1,1-^  ‘is 


(h)  ^(h) 


The  i*"^  element  of  the  vector  is  thus  the  contribution  to  the  (h)j 

row  and  (h)i^^  column  of  the  nvitrix.  The  integrals 


i “ / 

• J •'/  l_  \ 


may  be  found  by  using  the  local  coordinates, 


u = x-x 


V = y-y 


= ^*(h)l  ^ ^h)2  ^ ^h)3^^^ 
= ^y(h)i  ^ ^(h)2  >'(h)i)^^ 


* ■"  • '*1 


The  following  table  gives  the  most  nsetul  results  (10): 


i + j 

I . . 
11 

0 

A (area  of  triangle) 

1 

0 

2 

A / i i i ) . i i , 

--(UlVi  f U^V^,  + UjVj) 

3,4 

To  ^ 

5 

^ "2^'2  ^ 

where  u 

V 


i 

i 


’‘(h)i 

^h)i 


X 

c 

V 

' c 


and 

A = |det  i/2 
2 2 

If  k = Rq  e where  £(x,y)  is  a function  of  position,  we  mav  represent 
c(x,y)  by  piecewise  polynomial  functions,  and  the  integrals  can  be  found 
accordingly. 


( 3 ) Solving  the  linear  equations 

The  dimension  of  the  unknown  vector  u is  usually  very  large. 
Fortunately  the  matrix  G(I,J)  is  usually  sparse  and  banded.  Referring 
to  Fig.  13,  we  notice  that  the  nodes  (1,2, 3, 4)  are  only  connected  to 
the  Interior  nodes  (5,6, 7,8)  which  are  only  connected  to  the  interior 
(1,2, 3, 4)  and  (9,10,11,2)  etc.  Representing  these  nodal  values  of  the 
field  by  vectors, 
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(Uj.U^.U^.U^) 


''1  = 
V2  = 

V3  = 

V4  = 


we  can  represent  the  banded  matrix  by  the  vector  equations 

B.v.  , + A.v.  + C.v.  + d,  = 0 (25) 

1 1-1  11  1 i*-!  1 

The  matrix  equation  is  thus,  (referring  to  Fig.  13) 


GU 


0 

0 


(26) 


where  all  A^'s  are  A x A sparse  matrices,  C^'s  are  lower  triangular 
sparse  matrices  (frequently  diagonal)  and  B^'s  are  upper  triangular 
sparse  matrices  (frequently  diagonal);  = 0 and  the  vector  dj 

represent  the  contributions  from  the  boundary  values. 

The  solution  can  be  obtained  by  the  following  block-by-block 
elimination  method.  Assume 


Vi-1  ° ’^i''i  ®i  i = 1,2,.  . . 


(27) 


Eq.  (25)  can  be  rewritten  as 


- - (A^  + B^R.)"^C^v^^j  - (A.  + B^R^)"^  (B.;^  + d^)  (28) 
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I 

1 


fi 


a^B&tSBSSiiamsSiSSSSSI^lS^Stia: 


1 


'I 


Comparing  (28)  with  (27)  with  i = i+1,  we  get 


(29) 


and 


* »i'l>  <»1  = 1 * ■',> 


(30) 


The  above  two  recurrence  formulas  can  be  initiated  with  i = 1 in  Eq . (25), 
which  results  in 


''l 


(31) 


Comparing  (31)  with  Eq.  (27)  for  i = 2,  we  get 


R2  - -A^ 


-Sj  . -A^'dj 

Comparing  (32)  with  (29)  and  (30),  it  is  clear  that 


(32) 


Rl  = 0| 

, = 0 J 


(33) 


The  end  condition  at  i = N is 


Bv  ,+A.v  + d = 0 

N N-1  ^ N N 


(34) 


Substituting  Eq.  (27),  with  i = N,  into  (34)  we  get. 


N+1 


Ci'y) 


To  summarize,  we  have  the  following  procedure: 


(I)  Use  till'  Initial  conditions  (It) 

(II)  To  generate  K^,  s^  (i  = 1,7,..., N+1)  from  the  recurrence 
formulas  (29)  and  (30). 
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(Ill)  Use  the  end  condition  (35) 


(IV)  To  calculate  (i  = N,  N-1,  N-2,...,l)  from  (27). 

Applications 

(A)  Two  Dimensional  Scattering 

Equation  (29)  clearly  indicates  that  all  the  matrices  can  be 

computed  independent  of  the  boundary  values.  Therefore,  we  may  obtain 

a set  of  the  solutions  of  the  interior  problem  based  on  a set  of  specified 

boundary  values  from  the  same  R^'s.  The  vectors  s s are  dependent  on 

the  boundary  values  but  their  computations  are  not  as  time  consuming  as 

those  for  R.'s.  From  the  interior  solutions  the  normal  derivatives  for 
1 

each  Dirichlet  problem  can  be  found. 

Let  the  exterior  solution  be  expanded  in  a finite  series  of  cylindrical 
harmonics,  assuming  symmetry  with  respect  to  the  x-axis, 

u(r)  ~ 3 (hr)  ^ ^ ^ (36) 

n = 0 

and  the  interior  solution  be  expanded  into  N + 1 linearly  independent 
solutions  described  in  the  previous  section, 

N 

u(r)  = y'  b u (r)  r < a (37) 

^ - n n — 

n = 0 

where  a is  the  radius  of  the  circumscribing  circle,  and  the  normal  derivatives 
3U 

^ (a)  are  found  via  the  solutions  u (r).  To  find  the  coefficients 
dr  n 

a and  b , we  equate  (36)  and  (37),  and  their  normal  derivatives  at 
n n 

r = a,  to  get 
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.71/2  N 

I \ T'  ‘1  n>J  H (ka)}  cos  m((i  dO 

Jo  n=  0 " 

N ,7t/2 

= 52  ^ \ ^ (a,<{i)  cos  ra^)  d^»  (m  = 0,1,2,. ..,N) 

n = 0 " Jo  " 


(2) 

3H 


f 9U 

I (tt  (a,>)  + o cos  n$  - — (ka)}  cos  nut  d(J^ 

Jo  n = 0 " 

N tt/2  AU  (a,vM 

= ^ b f ■ - cos  m(l  di^  (m  = 0, 1 , . . . , N) 

n = 0 " Jo 


— ^ cos  m(l  di^  (m  = 0, 1 , . . . , N) 


The  above  should  result  in  a 2(N+1)  x 2(N+1)  matrix  from  which  the 

coefficients  can  be  found.  A word  of  advice  here  regarding  the  normal 
.,,.inc  ( 2 ) 


derivatives 


(r ) and  n 


(kr).  In  general,  the  exact  values  of  the 


normal  derivatives  of  the  incident  field  and  the  Hankel's  functions  are  known 
on  r = a,  however,  it  is  advisable  to  use  their  approximate  values 

— ^ and  — : (a)  instead.  The  reason  is  that  the  right  side  of  (39) 

Ar  Ar 

contains  appro.xinate  derivatives  obtained  from  the  finite  difference 
of  the  nodal  values.  They  are  not  equal  to  the  exact  derivatives.  The 
equality  of  the  approximate  and  the  exact  derivatives  can  be  enforced  in 
(39)  only  throufji the  alteration  of  the  coefficients  from  their  proper 
values.  Indeed,  equation  (39)  is  more  meaningful  if  both  sides  have  the 
same  order  of  approximation. 

A few  solutions  of  the  2-dimensional  scattering  Sfe  shown  in 
Fig.  15-16.  The  details  of  the  2-dimensional  solution  can  be  found  in 
Ref.  [3]. 


( B ) Three  Dimensiona 1 S cattering 
( 1 ) Formulation 

A truly  arbitrary  three  dimensional  scattering  problem  of  reasonable 
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size  is  not  going  to  be  solved  economically  in  the  foreseeable 
future  using  any  presently  available  technique.  The  three  dimensional 
scatterers  that  can  be  investigated  reasonably  economically  are  axially 
symmetrical.  An  axially  symmetrical  scatterer  with  non-axially  symmetric 
incident  fields  can  be  decomposed  into  azimuthal  modes. 


E(p,z,i(i)  = e^(o,z')e*'^’ 

m = 

OD 

riQH(p,z,4i)  = ^ ii^(p,z)e' 


where  Pq  = 120tt  ohms.  We  shall  use  the  cylind  cal  coordinates  for  the 
interior  problem  instead  of  spherical  coordinates  so  that  the  integrals 
resulting  from  the  finite  methods  will  be  more  manageable.  The  above 
azimuthal  components  can  be  solved  Independently  of  one  another. 

Let  R = kp  and  Z = kz,  and  p and  be  the  two  scalar  potentials, 


= 41  • e 


= 4,  . h (- 

R ^ m 

The  rest  of  the  field  components  can  be  derived  from  the  potentials 

9 

by  the  following  formulas: 


4i  X e = jf  (m4  x V4;  - Ru  V4; 

ra  m m,l  r m,2 


4 X ii  = jf  (mi  X Vij;  + Rc  V4  ) 
^ m m m,2  rm,l 


where 


f » Ip  (R,/.)i  (K,Z)R  - m‘| 
m r r 
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The  differential  equations  of  ij/  , and  ip  - are 

rn , 1 m , 2 


V • [f  (Rc  Vi|/  + raif  X V'^  1 + t +'  ,, 

m rm,l  m,2  rra.l 


/R  = 0 


V • [f  (Re  V’4;  - m;}>  X Vil  , ) 1 u 'I)  -/R  = 0 

m rra,2  ra,l  rm,2 


These  differential  equations  can  be  recast  into  the  minimization  of  the 


functional, 


F = J L(R,Z,'^^,'^,,VC,^,7i,,,)dRdZ 


L = f [Viii  • (Rc  Vij;  + mil  X Vij;  + Vi  (Ru^Vii)  . - mj  x 7v  ,) 
m m,l  r m,l  m,2  m,2  2 m,2  m,l 


r m, i r m, z 


It  is  noticed  that  f (R,Z)  is  singular  at  the  surfaces  where 
m 

r/u  e = |m|.  This  can  happen  when  p and  e are  both  real.  The  nodal 
r r ' ' r r 

fields  do  not,  however,  behave  erratically  at  these  surfaces  and,  in  fact, 
except  at  points  which  correspond  to  material  interfaces,  the  modal 
field  components  are  uniformly  holoraorphic. 


(2)  Singular  Integrals 

In  the  application  of  finite  metiiods  to  (47)  and  (48)  or  (50)  the 
singular  integrals  of  the  type. 


V,  ■ J 


R^Z" 

2 

2 e u R - m 
r r 


are  often  encountered.  The  integration  of  (51)  may  be  effected  as  follows: 


1 


-24- 


Consider  the  integral 


I = I g(R)Z  dRdZ  . 


Using  the  two  dimensional  Stoke 's  theorem. 


V X V • ds 


and  letting 


V(R,Z)  = - R . 


we  have 


V X V = -g(R)Z  4>  and  ds  = -<J)dRdZ 


and  (53)  becomes 


( g(R)Z®dRdZ  = g(R)Z®'^^  dR  . 

'n  ^ •'an 


Applying  (55)  to  (51)  we  get 


/ 2 2 
an  e U R - m 

r r 


For  a typical  element  the  clockwise  contour  of  (56)  is  shown  in  Fig.  17. 


Consider  a typical  integration  along  a line  segment,  such  as  from  the  i 


node  to  the  j node.  Along  this  path  the  linear  functional  dependence 


of  Z(R)  = oR  + B where 


(Z  - Z^) 

« = (Rj-'-R-.-)-  6 = Z^  - aR^ 


is  substituted  into  the  integrand  of  (56),  and  the  term  Z**"*"^  = (aR  + B)^^' 


is  expanded  into  binomial  series  in  powers  of  R.  The  resultant  integrals 


will  be  of  the  form. 


P 


f‘^  i r" 

'^2 

^ JR.  k“R  - ra" 

1 


whtre  n = 0,1,. ...r  + s + 1,  for  n ~ 0,1  wt‘  havo 


Ip(O.m)  = 


{ln(kR  - m)  - ln(kR  + m)  j ^ , n ^ 0 


1 (l,ra)  = - — , lln(kR  - m)  ln(kH  + tn)  ) • all  m (^0) 

^ 2k^  ^R. 

1 

Thfc  remainder  of  the  inteKr.ils  m.iy  be  generated  using  the  recur  renee  formula 


I (n,m)  = l^^(n,0)  + l|.^(n  - 2,m)  . (ol) 

< ^ k ' 

An  important  result  to  be  noted  Is  tlie  fact  Ifiat  even  for  the  ca^.e  of  lossless 

media,  where  k = p is  real,  the  quadrature  in  will  produce 

r r Q 

complex  values  wlien  the  ituegration  passes  through  the  simple  pole  singularity 
in  the  integrand.  The  pole  locations  for  the  case  of  real  and  complex  k 
are  illustrated  in  Fig.  18,  and  hence  the  integration  path  def  orir.at  ion 
necessary  for  real  k,  where  •'  -•  ().  in  pathlng  through  the  polo,  with 
increasing  R,  the  appropriate  natural  log  function  in  (.59)  and  (60) 
will  pick  up  an  additional  j :r  to  add  to  the  integral,  which  is  contributed 
by  the  residue  of  the  integral.  The  Global  System  matrix,  will  thus  include 
complex  elements  even  in  the  lossless  media  tase. 

Results 

Using  the  above  described  formulas  and  computational  techniques, 
we  are  able  to  compute  the  scattering  of  plane  waves  by  dielectric 
bodies  of  revolution  of  various  shapes. To  verify  the  computed  results 
we  first  compute  the  scattering  bv  a liielectric  sphere,  the  exact  result 


of  which  is  known-  in  t lie  computations,  we  deliberately  off  set  the  center 


L 


of  the  scatterer  from  tlie  origin,  so  tliat  computationally  speaking 


the  geometry  is  no  longer  regul.ir.  Mie  comp  ir Ison  of  the  results  Is 
shown  in  Fig.  19  for  tfie  scattering  amplitude  and  Fig.  20  for  the  phase. 

As  computational  techniques  becomes  mote  sophisticated,  verification 
of  the  computed  results  soon  becomes  a problem,  because  classical  solutions 
of  a few  special  geometries  will  not  be  sufficient  to  confirm  tlie  comput  .at  Ions 
of  complex  problems.  ExpcrininUal  vt'r  it  ication  may  eventually  be  tlie  only 
way  to  confirm  the  validity  of  general  programs.  We  have  both  computed 
and  measured  the  bistatic  scattering  of  a finite  dielectric  cylinder  with 
a spherical  void  in  the  center  such  as  shown  in  Fig.  21.  The  computed 
and  measured  results  are  shown  in  Fig.  22  and  21.  Their  agreement  is 
truly  remarkable.  These  results  should  definitively  confirm  the  applicability 
of  the  unimoment  method  in  electromagnetic  scattering  problems. 


Further  development 

At  the  present  using  line.ir  finite  elements  methods,  it  Is  possible 
to  solve  axially  s>'nnnet  r ica  I dielectric  srarieTing  problems  wiiich  ti.ave  a 
maximum  dimension  of  about  using  the  tdJC  7b00  computer.  Actually 

it  is  difficult  to  he  exact  regarding  tile  limitation  of  the  method  because 
it  should  depend  on  tiie  dl elect  tic  constant  and  geometry  of  the  scatterer. 
At  this  stage,  however,  it  is  quite  evident  that  tiie  Immediate  limitation 
of  the  method  is  the  capacity  ot  tiie  computer  in  handling  finite  element 

j equations  in  a large  closed  region.  Fortunately,  finite  element  tecimiques 

i 

j have  been  extensively  investigated  by  st rue tural  engineers  and  improvements 

I of  these  techniques  have  been  reported  at  each  of  the  related  conferences. 

I 

I The  finite  element  teciinique  desrrihed  in  this  paper  is  quite  elementary. 

I Yet,  with  this  elementarv  appro. uti,  we  have  .already  extended  the  art  of 


I 
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scattering  computation  quite  significantly.  It  is  conceivable  tliat  more 
sophisticated  finite  metliods  could  enh.ince  tlie  dimensions  of  solvable 
interior  problems  to  the  point  tli.it  tiie  size  of  tlie  matrix,  resulting  from 
enforcing  the  continuity  conditions  on  tlie  separable  surface,  becomes  the 
limiting  factor  oi  the  computation.  Therefore,  the  most  needed  effort 
is  to  exploit  the  power  of  the  finite  element  method  and  to  popularize  Its 
use  .iraong  the  electromagnetic  community. 

The  application  of  the  technique  is  very  broad  indeed.  Current 
problems,  such  as,  scattering  by  advanced  composite,  m.iterlals  absorption  by 
biological  material,  scattering  by  buried  obstacles,  just  to  name  a few, 
can  all  be  analyzed  with  finite  metliods.  It  is  also  conceivable  that 

the  unimoment  method  can  he  ccimbined  with  the  method  of  moments  to  solve 
a variety  of  problems,  which  consist  of  wires  and  hodles  of  revolution. 

In  conclusion,  we  would  like  to  (mphasize  tint  the  .•<ppl  ieat  ion  of  finite 
methods  in  antenna  and  scattering  problems  is  not  an  alternative  to  the 
method  of  moment,  rather  it  is  its  supplement.  It  is  obvious  that  the 
finite  method  is  most  uneconomical  where  the  method  of  moments  is  most 
economical, such  as  in  thin  wire  problems,  wliereas  the  finite  method  is 
most  attractive  where  the  method  of  moments  is  most  clumsy, such  as  scattering 
by  inhomogeneous  material  bodies. 
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Fig.  11.  (a)  Left  handed  and  (b)  right  handed  elements. 


Fig.  14.  A linear  weighting  function. 
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Fig.  18.  Tlie  pole  iocat.ions  of  real  .iiul  complex  k. 


Fig.  20.  The  comparison  of  phases  of  scattered  fields  by  numerical  and 
analvtica]  solutions  of  a dielectric  spiiere. 


Fig.  22.  Computed  and  measured  amplitudes  of  scattc’red  fields  of  a 

hollow  dielectric  cylinder  illuminated  by  a 90°  incident  plane 
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